#rm(list=ls())
#setwd("/Main Text-Figures/Figure3/")

library(pscl)

########################################

alpha <- .5
k <- 0.2
N <- 100
v <- 500
nsim <- 500

#heterogeneity scores
d1 <- seq(0.5,2, .1)
d2 <- d1

cohesion.right.all.rprop.sim <- NULL
cohesion.right.rcv.rprop.sim <- NULL
cohesion.right.all.lprop.sim <- NULL
cohesion.right.rcv.lprop.sim <- NULL

cohesion.left.all.rprop.sim <- NULL
cohesion.left.rcv.rprop.sim <- NULL
cohesion.left.all.lprop.sim <- NULL
cohesion.left.rcv.lprop.sim <- NULL


cohesion.right.all.sim <- NULL
cohesion.left.all.sim <- NULL
cohesion.right.rcv.sim <- NULL
cohesion.left.rcv.sim <- NULL
diff.left.cohesion.sim <- NULL
diff.right.cohesion.sim <- NULL

for(z in 1:nsim){
  phi <- runif(1,.2,.8)
  
  b1 <- runif(v,0,1)
  b2 <- b1
  sq1 <- runif(v,0,1)
  sq2 <- sq1
  
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  
  lco.lprop.all <- NULL
  lco.rprop.all <- NULL
  lco.lprop.rcv <- NULL
  lco.rprop.rcv <- NULL
  
  rco.lprop.all <- NULL
  rco.rprop.all <- NULL
  rco.lprop.rcv <- NULL
  rco.rprop.rcv <- NULL
  
  diff.right.cohesion <- NULL
  diff.left.cohesion <- NULL
  cohesion.left.all <- NULL
  cohesion.left.rcv <- NULL
  cohesion.right.all <- NULL
  cohesion.right.rcv <- NULL
  
  
  for(J in 1:length(d1)){
    
########################################
    
    hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1[J],length(b1)),b1,sq1)
    colnames(hypo.pro.1d) <- c("dim","d","b","sq")
    hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2[J],length(b2)),b2,sq2)
    colnames(hypo.pro.2d) <- c("dim","d","b","sq")
    hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
    hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
    
########################################
    
    xL1 <- c(0,runif(npL-1,min=0-d1[J],max=d1[J])) #Going to use this that i've put the party leaders in the first and last spots
    xL2 <- xL1
    xR1 <- c(runif(npR-1,min=1-d1[J],max=1+d1[J]),1) #Going to use this that i've put the party leaders in the first and last spots
    xR2 <- xR1
    
    left.party <- as.data.frame(cbind(xL1,xL2,rep(0,npL)))
    colnames(left.party) <- c("x1","x2","partyid")
    right.party <- as.data.frame(cbind(xR1,xR2,rep(1,npR)))
    colnames(right.party) <- c("x1","x2","partyid")
    all.leg <- as.data.frame(rbind(left.party,right.party))
    
########################################
#who votes for it?

    vote.func1 <- function(x1,b1,sq1){
      as.numeric(abs(x1-b1)  < abs(x1-sq1) )
    }
    vote.func2 <- function(x2,b2,sq2){
      as.numeric(abs(x2-b2)  < abs(x2-sq2) )
    }

    out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
    for(i in 1:N){
      for(j in 1:v){
        out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
      }
    }
    for(i in 1:N){
      for(j in (v+1):(2*v)){
        out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
      }
    }
    
########################################
#does it pass?
    
    passes <- as.numeric(colSums(out) > .5*N)
    hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
    
    for (i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
    }
########################################
#unity scores?
    
    left.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
    }
    right.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
    }

    out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] - right.aligned[j]
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]-left.aligned[j]
        } else {
          0
        }
      }
    }
    
    out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] 
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]
        } else {
          0
        }
      }
    }
########################################
#rcv requested?
    
    out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
    for(i in 1:N){
      for(j in 1:(v*2)){
        out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
      }
    }
    hypo.pro.all$left.cscore <- out.cohesion[1,] 
    hypo.pro.all$right.cscore <- out.cohesion[N,]
    hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
    hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
    hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
    hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
    
    hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)

########################################
#would this ever get proposed?
    
    for(i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$proposing.party[i] <- if(hypo.pro.all$b[i] < hypo.pro.all$sq[i]) 0 else 1 
    }
    
    hypo.pro.all$proposed.left <- rep(0,length(nrow(hypo.pro.all)))
    for(i in 1:nrow(hypo.pro.all)){
      if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else {
        0
      }
    }
    hypo.pro.all$proposed.right <- rep(0,length(nrow(hypo.pro.all)))
    for(i in 1:nrow(hypo.pro.all)){
      if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else {
        0
      }
    }
  
    hypo.pro.all$proposed <- hypo.pro.all$proposed.left+hypo.pro.all$proposed.right

    all.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1),]
    all.proposals <- all.proposals[which(all.proposals$rcv.request==0),]
      all.proposals.left <- all.proposals[which(all.proposals$proposed.left==1),]
      all.proposals.right <- all.proposals[which(all.proposals$proposed.right==1),]
    rcv.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1 & hypo.pro.all$rcv.request==1),]
      rcv.proposals.left <- rcv.proposals[which(rcv.proposals$proposed.left==1),]
      rcv.proposals.right <- rcv.proposals[which(rcv.proposals$proposed.right==1),]
########################################

    vote.matrix.all <- out[,all.proposals$id]
      vote.matrix.all.left <- out[,all.proposals.left$id]
      vote.matrix.all.right <- out[,all.proposals.right$id]
    vote.matrix.rcv <- out[,rcv.proposals$id]
      vote.matrix.rcv.left <- out[,rcv.proposals.left$id]
      vote.matrix.rcv.right <- out[,rcv.proposals.right$id]
    
########################################
#Cohesion
      
    cpd.all <- as.data.frame(c(all.proposals$left.cscore.abs,all.proposals$right.cscore.abs))
    cpd.all$id <- c(rep("L",length(all.proposals$left.cscore.abs)),rep("R",length(all.proposals$right.cscore.abs)))
    colnames(cpd.all) <- c("cscore","id")
    
    cpd.rcv <- as.data.frame(c(rcv.proposals$left.cscore.abs,rcv.proposals$right.cscore.abs))
    cpd.rcv$id <- c(rep("L",length(rcv.proposals$left.cscore.abs)),rep("R",length(rcv.proposals$right.cscore.abs)))
    colnames(cpd.rcv) <- c("cscore","id")
    
    cpd.right <- as.data.frame(c(all.proposals$right.cscore.abs,rcv.proposals$right.cscore.abs))
    cpd.right$id <- c(rep("ALL",length(all.proposals$right.cscore.abs)),rep("RCV",length(rcv.proposals$right.cscore.abs)))
    colnames(cpd.right) <- c("cscore","id")

    cpd.left <- as.data.frame(c(all.proposals$left.cscore.abs,rcv.proposals$left.cscore.abs))
    cpd.left$id <- c(rep("ALL",length(all.proposals$left.cscore.abs)),rep("RCV",length(rcv.proposals$left.cscore.abs)))
    colnames(cpd.left) <- c("cscore","id")

    cohesion.left.temp <- mean(all.proposals$left.cscore.abs)
    cohesion.left.rcv.temp <- mean(rcv.proposals$left.cscore.abs)
    cohesion.left.all <- c(cohesion.left.all,cohesion.left.temp)
    cohesion.left.rcv <- c(cohesion.left.rcv,cohesion.left.rcv.temp)
    
    cohesion.right.temp <- mean(all.proposals$right.cscore.abs)
    cohesion.right.rcv.temp <- mean(rcv.proposals$right.cscore.abs)
    cohesion.right.all <- c(cohesion.right.all,cohesion.right.temp)
    cohesion.right.rcv <- c(cohesion.right.rcv,cohesion.right.rcv.temp)
    
    diff.L <- mean(all.proposals$left.cscore.abs) - mean(rcv.proposals$left.cscore.abs)
    diff.left.cohesion <- c(diff.left.cohesion,diff.L)
    diff.R <- mean(all.proposals$right.cscore.abs) - mean(rcv.proposals$right.cscore.abs)
    diff.right.cohesion <- c(diff.right.cohesion,diff.R)
    
    lco.lprop.all.temp <- mean(all.proposals.left$left.cscore.abs)
    lco.rprop.all.temp <- mean(all.proposals.right$left.cscore.abs)
    lco.lprop.all <- c(lco.lprop.all,lco.lprop.all.temp)
    lco.rprop.all <- c(lco.rprop.all,lco.rprop.all.temp)
    
    lco.lprop.rcv.temp <- mean(rcv.proposals.left$left.cscore.abs)
    lco.rprop.rcv.temp <- mean(rcv.proposals.right$left.cscore.abs)
    lco.lprop.rcv <-  c(lco.lprop.rcv,lco.lprop.rcv.temp)
    lco.rprop.rcv <- c(lco.rprop.rcv,lco.rprop.rcv.temp)

    rco.lprop.all.temp <- mean(all.proposals.left$right.cscore.abs)
    rco.rprop.all.temp <- mean(all.proposals.right$right.cscore.abs)
    rco.lprop.all <- c(rco.lprop.all,rco.lprop.all.temp)
    rco.rprop.all <- c(rco.rprop.all,rco.rprop.all.temp)
    
    rco.lprop.rcv.temp <- mean(rcv.proposals.left$right.cscore.abs)
    rco.rprop.rcv.temp <-  mean(rcv.proposals.right$right.cscore.abs)
    rco.lprop.rcv <- c(rco.lprop.rcv,rco.lprop.rcv.temp)
    rco.rprop.rcv <- c(rco.rprop.rcv,rco.rprop.rcv.temp)    
  }
  
  diff.left.cohesion.sim <- rbind(diff.left.cohesion.sim,diff.left.cohesion)
  diff.right.cohesion.sim <- rbind(diff.right.cohesion.sim,diff.right.cohesion)

  cohesion.right.all.sim <- rbind(cohesion.right.all.sim,cohesion.right.all)
  cohesion.right.rcv.sim <- rbind(cohesion.right.rcv.sim,cohesion.right.rcv)
  cohesion.left.all.sim <- rbind(cohesion.left.all.sim,cohesion.left.all)
  cohesion.left.rcv.sim <- rbind(cohesion.left.rcv.sim,cohesion.left.rcv)
  
  cohesion.right.all.rprop.sim <- rbind(cohesion.right.all.rprop.sim, rco.rprop.all)
  cohesion.right.all.lprop.sim <- rbind(cohesion.right.all.lprop.sim, rco.lprop.all)
  cohesion.right.rcv.rprop.sim <- rbind(cohesion.right.rcv.rprop.sim, rco.rprop.rcv)
  cohesion.right.rcv.lprop.sim <- rbind(cohesion.right.rcv.lprop.sim, rco.lprop.rcv)
  
  cohesion.left.all.lprop.sim <- rbind(cohesion.left.all.lprop.sim,lco.lprop.all)
  cohesion.left.all.rprop.sim <- rbind(cohesion.left.all.rprop.sim,lco.rprop.all)
  cohesion.left.rcv.lprop.sim <- rbind(cohesion.left.rcv.lprop.sim,lco.lprop.rcv)
  cohesion.left.rcv.rprop.sim <- rbind(cohesion.left.rcv.rprop.sim,lco.rprop.rcv)
}

rcohesion.rprop.all <- NULL
for(i in 1:ncol(cohesion.right.all.rprop.sim)){
  rcohesion.rprop.all[i] <- mean(cohesion.right.all.rprop.sim[,i],na.rm=TRUE)
}

rcohesion.rprop.rcv <- NULL
for(i in 1:ncol(cohesion.right.rcv.rprop.sim)){
  rcohesion.rprop.rcv[i] <- mean(cohesion.right.rcv.rprop.sim[,i],na.rm=TRUE)
}

rcohesion.lprop.all <- NULL
for(i in 1:ncol(cohesion.right.all.lprop.sim)){
  rcohesion.lprop.all[i] <- mean(cohesion.right.all.lprop.sim[,i],na.rm=TRUE)
}

rcohesion.lprop.rcv <- NULL
for(i in 1:ncol(cohesion.right.rcv.lprop.sim)){
  rcohesion.lprop.rcv[i] <- mean(cohesion.right.rcv.lprop.sim[,i],na.rm=TRUE)
}

lcohesion.rprop.all <- NULL
for(i in 1:ncol(cohesion.left.all.rprop.sim)){
  lcohesion.rprop.all[i] <- mean(cohesion.left.all.rprop.sim[,i],na.rm=TRUE)
}

lcohesion.rprop.rcv <- NULL
for(i in 1:ncol(cohesion.left.rcv.rprop.sim)){
  lcohesion.rprop.rcv[i] <- mean(cohesion.left.rcv.rprop.sim[,i],na.rm=TRUE)
}

lcohesion.lprop.all <- NULL
for(i in 1:ncol(cohesion.left.all.lprop.sim)){
  lcohesion.lprop.all[i] <- mean(cohesion.left.all.lprop.sim[,i],na.rm=TRUE)
}

lcohesion.lprop.rcv <- NULL
for(i in 1:ncol(cohesion.left.rcv.lprop.sim)){
  lcohesion.lprop.rcv[i] <- mean(cohesion.left.rcv.lprop.sim[,i],na.rm=TRUE)
}

mean.cohesion.right.all <- NULL
for(i in 1:ncol(cohesion.right.all.sim)){
  mean.cohesion.right.all[i] <- mean(cohesion.right.all.sim[,i],na.rm=TRUE)
}

mean.cohesion.right.rcv <- NULL
for(i in 1:ncol(cohesion.right.rcv.sim)){
  mean.cohesion.right.rcv[i] <- mean(cohesion.right.rcv.sim[,i],na.rm=TRUE)
}

mean.cohesion.left.all <- NULL
for(i in 1:ncol(cohesion.left.all.sim)){
  mean.cohesion.left.all[i] <- mean(cohesion.left.all.sim[,i],na.rm=TRUE)
}
mean.cohesion.left.rcv <- NULL
for(i in 1:ncol(cohesion.left.rcv.sim)){
  mean.cohesion.left.rcv[i] <- mean(cohesion.left.rcv.sim[,i],na.rm=TRUE)
}

########################################

load("cohesion-fig3ab-data.RData")

#pdf("Fig3a_.pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,3,1))
plot(x=d1,y=lcohesion.lprop.all,type='l',main="",xlab="",ylab="Party Cohesion",yaxt='n',
     lwd=2,col="black",ylim=c(.5,1))
lines(x=d1,y=lcohesion.rprop.all,lty=1,col="grey60",lwd=2)
lines(x=d1,y=lcohesion.lprop.rcv,col="black",lty=3,lwd=2)
lines(x=d1,y=lcohesion.rprop.rcv,lty=3,col="grey60",lwd=2)
axis(side=2,at=seq(.3,1,.1),labels=seq(.3,1,.1),las=1)
mtext(side=3,line=1,"(a) Left Party",cex=1.25)
mtext(side=1,line=2.5,"d: Party Heterogeneity")

par(mar=c(4,1,3,4))
plot(x=d1,y=rcohesion.rprop.all,type='l',ylim=c(.5,1),main="",xlab="",ylab="",yaxt='n',col="grey60",lwd=2)
lines(x=d1,y=rcohesion.rprop.rcv,lty=3,col="grey60",lwd=2)
lines(x=d1,y=rcohesion.lprop.all,col="black",lwd=2)
lines(x=d1,y=rcohesion.lprop.rcv,lty=3,col="black",lwd=2)
axis(side=4,at=seq(.3,1,.1),labels=seq(.3,1,.1),las=1)
mtext(side=3,line=1,"(b) Right Party",cex=1.25)
mtext(side=1,line=2.5,"d: Party Heterogeneity")
mtext(side=4,"Party Cohesion",line=3)
#dev.off()
